%Test the simiImplicitmove for moving mesh
clear all;
close all;

xl = 1;
yl = 1;
global nelx;
nelx = 64;
nely = 64;
dkesi = xl/nelx;
dyita = yl/nely;
[kesi,yita] = meshgrid(1:nelx, 1:nely);
kesi = kesi/nelx;
yita = yita/nely;
x = sin(pi/2*kesi);
y = yita;    % transform the mesh
s = zeros(nelx,nely);
%figure(1), surf(kesi,yita,s); % computational mesh
%figure(2), surf(x,y,s); % test the mesh

ne = nelx * nely;
kappa = 5;
pdecycle = 50;
dt = .5*dkesi;
phi = ones(nelx,nely);
R = .3; % radius for initial circle

for j = 1 : nely
    for i = 1 : nelx
    if (x(i,j)-xl/2)^2+(y(i,j)-yl/2)^2 < R^2
        phi(i,j) = 0;
    end
   end
end
figure(3), surf(x,y,phi); % initial shape
axis equal; axis tight; axis off;pause(1);%Above 3 lines are plotting phi. 

phi=reshape(phi, nelx*nely, 1);
mphi = phi - phi.^3; 

xkesi = zeros(size(x));
ykesi = zeros(size(x));
xyita = zeros(size(x));
yyita = zeros(size(x)); % solve for x_yita, x_kesi, y_yita, y_kesi

for i = 2 : nelx-1
    for j = 2 : nely-1
    xkesi(i,j) = (x(i+1,j)-x(i-1,j))/(2*dkesi);
    xyita(i,j) = (x(i,j+1)-x(i,j-1))/(2*dyita);
    ykesi(i,j) = (y(i+1,j)-y(i-1,j))/(2*dkesi);
    yyita(i,j) = (y(i,j+1)-y(i,j-1))/(2*dyita);
    end
    xkesi(i,1) = (x(i+1,1)-x(i-1,1))/(2*dkesi);
    xyita(i,1) = (x(i,2)-x(i,1))/dyita;
    ykesi(i,1) = (y(i+1,1)-y(i-1,1))/(2*dkesi);
    yyita(i,1) = (y(i,2)-y(i,1))/dyita;
    
    xkesi(i,nely) = (x(i+1,nely)-x(i-1,nely))/(2*dkesi);
    xyita(i,nely) = (x(i,nely)-x(i,nely-1))/dyita;
    ykesi(i,nely) = (y(i+1,nely)-y(i-1,nely))/(2*dkesi);
    yyita(i,nely) = (y(i,nely)-y(i,nely-1))/dyita;
end

for j = 2 : nely-1
    xkesi(1,j) = (x(2,j)-x(1,j))/dkesi;
    xyita(1,j) = (x(1,j+1)-x(1,j-1))/(2*dyita);
    ykesi(1,j) = (y(2,j)-y(1,j))/dkesi;
    yyita(1,j) = (y(1,j+1)-y(1,j-1))/(2*dyita);
    
    xkesi(nelx,j) = (x(nelx,j)-x(nelx-1,j))/dkesi;
    xyita(nelx,j) = (x(nelx,j+1)-x(nelx,j-1))/(2*dyita);
    ykesi(nelx,j) = (y(nelx,j)-y(nelx-1,j))/dkesi;
    yyita(nelx,j) = (y(nelx,j+1)-y(nelx,j-1))/(2*dyita);
end
    xkesi(1,1) = (x(2,1)-x(1,1))/dkesi;
    xyita(1,1) = (x(1,2)-x(1,1))/dyita;
    ykesi(1,1) = (y(2,1)-y(1,1))/dkesi;
    yyita(1,1) = (y(1,2)-y(1,1))/dyita;
    
    xkesi(1,nely) = (x(2,nely)-x(1,nely))/dkesi;
    xyita(1,nely) = (x(1,nely)-x(1,nely-1))/dyita;
    ykesi(1,nely) = (y(2,nely)-y(1,nely))/dkesi;
    yyita(1,nely) = (y(1,nely)-y(1,nely-1))/dyita;
    
    xkesi(nelx,1) = (x(nelx,1)-x(nelx-1,1))/dkesi;
    xyita(nelx,1) = (x(nelx,2)-x(nelx,1))/dyita;
    ykesi(nelx,1) = (y(nelx,2)-y(nelx,1))/dkesi;
    yyita(nelx,1) = (y(nelx,2)-y(nelx,1))/dyita;
    
    xkesi(nelx,nely) = (x(nelx,nely)-x(nelx-1,nely))/dkesi;
    xyita(nelx,nely) = (x(nelx,nely)-x(nelx,nely-1))/dyita;
    ykesi(nelx,nely) = (y(nelx,nely)-y(nelx-1,nely))/dkesi;
    yyita(nelx,nely) = (y(nelx,nely)-y(nelx,nely-1))/dyita;
 
    J = xkesi*yyita - xyita*ykesi;   %approximate for J
    
for i = 1 : pdecycle 
%[phi,A,p,q,r] = simiImplicitmove(phi,mphi,nelx,nely,xkesi,ykesi,xyita,yyita,J,dkesi,dyita,dt,kappa);
[phi] = simiImplicitmove2(phi,mphi,nelx,nely,xkesi,ykesi,xyita,yyita,J,dkesi,dyita,dt,kappa);
phi =  max(0,min(1, phi));

figure(4), surf(x,y,1-phi),axis([0 1 0 1]), pause(.1); 
%plotting phi. 

phi = reshape(phi,ne,1);
mphi = phi - phi.^3;
end    

 






